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Abstract 

We relate the quantum dynamics of the Bose-Hubbard model (BHM) to the 
semiclassical nonlinear equations that describe an array of interacting Bose conden- 
sates by implementing a standard variational procedure based on the coherent state 
method. We investigate the dynamics of the two-site BHM from the purely quan- 
tum viewpoint by recasting first the model within a spin picture and using then the 
related dynamical algebra. The latter allows us to study thoroughly the energy spec- 
trum structure and to interpret quantally the classical symmetries of the two-site 
dynamics. The energy spectrum is also evaluated through various approximations 
relying on the coherent state approach. 

1 Introduction 

The simplest way to describe an array of interacting bosonic wells is represented by 
the second quantized Hamiltonian for charged bosons 

T 

H bh = J2 \-~^ n i + Uni( ni (afaj + a+a^j , (1) 

* (iJ) 

where the operators = afa>i count the number of bosons at the lattice site i, 
and the destruction (creation) operators ai, (af) obey the canonical commutation 
relations [a*, aj] = The parameters T, U > 0, and \i represent the hopping am- 
plitude, the strength of the Coulomb on-site repulsion, and the chemical potential, 
respectively. 



Hamiltonian (|l]), known in the literature Q as the Bose-Hubbard Model (BHM), 
was originally derived as the lattice version of the second quantized |^| 4 -field theory 
used to model superfluids adsorbed in porous media. It describes in fact a gas of 
identical charges that hop through a D dimensional lattice, and exhibit a repulsive 
action whenever two (or more) charges are situated at the same site. The bosonic 
character of the lattice stands in the fact that, in principle, an unlimited number 
of charges is allowed to occupy each lattice sites. On this account the BHM is able 
to mimic appropriately various features of systems such as granular superconduc- 
tors, short-length superconductors and Josephson junction arrays, hence deserving 
a growing attention in the recent years (see references in Ref. (2|). For the same 
reason the BHM seems to supply the appropriate theoretic framework in which 
describing an array of interacting Bose-Einstein condensates. 

The recent years have been characterized by remarkable progresses in construct- 
ing experimental devices able to confine some hundreds of bosons so as to form real 
Bose-Einstein condensates (BEC) ||. This has prompted as well the investigation 
of more structured situations in which two condensates (regarded as sites or wells 
filled by condensed bosons within the BHM picture) can interact via tunneling pro- 
cesses 0. In particular, a quite rich scenery of dynamical behaviors has been found 
to characterize the two-site case (named boson Josephson junction) in Ref. ||, 
where two coupled Gross-Pitaevskii equations for two condensate wavefunctions are 
assumed as the minimal interaction scheme whereby describing interwell processes. 
This represents however the most general array of coupled BEC's so far realized 
even if BEC's with complex architectures have been constructed 0, and a large 
amount of both experimental and theoretic work has been devoted to construct ar- 
rays within the physics of uncondensed ultracold atoms in optical lattices (see Refs. 
0, (8| and references therein). 

Despite its almost purely theoretic value, it is interesting to illustrate for con- 
ceptual reasons the link between the BHM and a system of many bosonic wells 
regarded as interacting Bose condensates. We reconstruct such a correspondence in 
Sec. 2 for the general case and employ it in Sec. 3 to show how the two-site BHM 
coincides, in a semiclassical picture, with a two-condensate system. 

In this paper we shall focus our attention on the quantum dynamics of the two- 
site case. Classically, the problem is known to be integrable since the dynamics 
of the complex site variables z±, z% required to account for the system state (these 
identify with the wavefunctions of the two-condensate system) can be recast in terms 
of two action-angle variables by exploiting the constant of motion M = \z\\ 2 + 
\z2 1 2 - The reliability of the classical picture comes, of course, from the fact that 
indeed the operators ng counts a large number of particle and thus identifies with 
its semiclassical counterpart |z^| 2 via the correspondence scheme discussed in Sec. 
2. 
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A smaller number of particles per site makes reasonable recovering the quantum 
picture both because this is experimentally possible, and because it is interesting 
to investigate thoroughly the changes in the structure of the quantum states when 
going from the pure quantum regime to situations with an intermediate, more or 
less pronouced, semiclassical character. Moreover, it is important to recall that 
the BHM exhibits a zero temperature [i-T phase diagram which contains regions 
of finite size (lobes) where the system loses its superfluid character becoming an 
insulator || . In this respect we are interested in making evident the cause of the 
presence of insulator regions in the phase diagram in relation with the structure of 
the two-site model eigenstates. 

In Sec. 3 we review the semiclassical form of the two-site dynamics and im- 
plement the standard requantization procedure for comparing it to the spectrum 
provided by the exact quantum treatment. In order to study the quantum dy- 
namics of the two-site BHM in Sec. 4 we represent the Hamiltonian within a spin 
algebra realization based on the generators of su(2), in which the index of the al- 
gebra representation is proportional to the total boson number N. This establishes 
the equivalence between the usual description of the BHM Hilbert space through 
the bosonic Fock basis and a simplified picture relying on the angular momentum 
eigenstates. The latter, in turn, allows one to use the su(2) coherent states (where 
the conservation of the total boson number has been already incorporated) when 
performing the classical limit. 

The spin picture also provides the basis for developing in Sec. 5 the single-boson 
picture of the system dynamics. This amounts to reconstructing the exact dynami- 
cal algebra |J of the two-site BHM. Explicitly, the model is reformulated within an 
enlarged algebra -which identifies with su(N+l) in its fundamental representation- 
where the two-site Hamiltonian takes the form of a many-boson Hamiltonian on a 
(N+l)-site lattice whose population is formed by a unique (effective) boson. The 
algebraic formulation we introduce is viable to BHM's with more than two site 
and provides a deeper insight of the dynamics complexity from the group-theoretic 
standpoint. Moreover, it appears particularly advantageous in studying the quan- 
tum counterpart of the nonintegrable classic dynamics of the N-site BHM (N > 2). 
This problem will be considered elsewhere. 

At the operational level, such an approach entails a reformulation of the quan- 
tum dynamics in terms of a (classical) Hamiltonian which exhibits a quadratic 
dependence on the components of the su(2) picture states. In particular, such an 
Hamiltonian is employed in Sec. 6 to analyse various dynamical regimes and to rec- 
ognize the distinctive characters of both the superfluid and the insulator behavior. 
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2 Arrays of condensates within the BHM pic- 
ture 

The link that relates the BHM to a system of interacting boson condensates is easily 
established by implementing the method used in Refs. [l(J for investigating the 
quantum dynamics of many-particle lattice models. To this end the one-dimensional 
(ID) Hamiltonian derived from (|l]) 



3=1 



H c = ^2 -^ n j + Unj(nj - 1) - — (a+a j+1 + a+ +1 cij 



(2) 



suffices to generate a chain of M interacting wells namely the generalization of the 
case with two interacting condensates. 

The method just mentioned combines the generalized coherent state picture of 
quantum systems with the time-dependent variational principle (TDVP) B, and 
allows one to reformulate the quantum dynamical problem relative to (g) by a 
semiclassical hamiltonian dynamics in terms of the expectation values 

zj := (Z\ aj \Z) , z* := {Z\a}\Z) (3) 

and generated within the TDVP procedure by Hamiltonian 7i c (Z, Z*) = (Z\H C \Z). 
The so-called trial macroscopic wave function \Z) here is assumed to have the form 
\Z) = Ylj \zj), where 1 < j < M and \zj) are the Glauber coherent states Q solu- 
tions of the standard equation aj\zj) = Zj \zj) for each j. The dynamical equations 
for Zj (those for z* are obviously obtained by complex conjugation) 

T 

ihzj = (2U\zj\ 2 - n)zj - - z i ( 4 ) 

are derived [§] from the semiclassical Hamiltonian 

T 



H C (Z, Z*) = 



E 

3 



V + U\z 3 \ 2 )\ Zj \ 2 - - (z*z j+1 + z* 



2 



(5) 



through the standard formulas Zj = {zj,Ti. c } based on the canonical Poisson brackets 
{z*,zc} = (i/%)6j t £. It is easily checked as well that the basic feature [H c , N] = 
of the quantum system has a classical counterpart given by {TC c ,j\f} = 0, where 
j\f = J2e \ z i\ 2 - The identification of Zj with the wavefunction of the j-th condensate 
establishes the link claimed above, showing explicitly that the effective dynamics of 
BHM generates the classical dynamics of coupled condensates. 
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3 Semiclassical picture of the two- well system 



It is interesting to derive the two-well dynamics from the general case of the M- 
site lattice described by Hamiltonian (|5|). Upon assuming z 2s = z 2 and z 2s +i = 
z\ (suppose that M is even and the lattice is endowed with periodic boundary 
condition), (||) leads to the two-site Hamiltonian 

H:=H c /(M/2) = \u{\z 1 \ A + \z 2 \ A )-li{\z 1 \ 2 + \z 2 \ 2 )-T{zlz2 + z* 2 z 1 )\ , (6) 

while (H) reduce to the two equations 

%hz\ = (2U\zi\ 2 — [i)Z\ — Tz 2 

(7) 

ihz 2 = (2U\z 2 \ 2 - n)z 2 - Tz\ 

This amounts to considering a special solution characterized by two independent 
complex variables with |zi| 2 + |z2| 2 representing the constant of motion. Remarkably, 
the same equations are obtained independently when assuming Ti to be the two-site 
Hamiltonian equipped with the Poisson brackets {z a ,z%} = 5 a b/ih, a,b = 1,2. It 
is easily recognized that g-site models, where qp = M for some p £ N, can be 
achieved by following the same procedure (z a+ k q = £ Q , < /c < p, ^ < a < q) 
and that their dynamics can always interpreted as solutions of the M-site model 
exhibiting the discrete translation symmetry Zj+^g = z j > Vj, < k < p. The 
physical interpretation of such a symmetry is given by the Fourier modes description 
of the dynamical variables \f~Mb q = T,jZjexp[iqj], where q := 2-Kq/M, j, q G [1, M], 
and M5(j — £) = H q exp[iq(j — £)]. By inserting the condition z 2 = Z2 S and z\ = Z2s+i 
one finds 6q = yM(zi + z 2 )/2, and bM/2 = \fM(z\ — z 2 )/2 that imply 



zi = (b + b M/2 )/(2VM) , z 2 = (b - b M /2)/(2VM) . 

Therefore the two-site model represents the simplest possible way to construct an 
excited state with respect to the condensed macroscopic state characterizing the 
zero temperature regime where Zj = bo/\^M, Vj. Hamiltonian (0) can be rewritten 
as 

'U. r0 U 



H = 



2 -M 2 - fiM + -D 2 - TVM 2 - D 2 cos(26) 



(8) 



where we have introduced the new canonical variables J\f = \z\\ 2 + \z2\ 2 , D = 
\zi\ 2 — \z2\ 2 , 6 = (0i — 02)/2, and ip = ((pi + <f>2)/2 implied by Zj = \zj\exp(i(j)j). 
The new coordinates are equipped with the Poisson brackets 

{9,M} = = {ip, D} , {6, D} = -1/h = {ip,N} , 
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that provide the equation of motions 



j 1 n 

D = 2TVN 2 -D 2 sin(20) , 6 = -UP- , = cos(29) , (9) 

VTV 2 - D 2 

the dot denotes the time derivative d/dr with r := t/Ti. Based on the complete 
integrability of the system, and proceeding along the same lines of Ref. |1C], one 
can reduce the above two equations to a unique one 

D 2 = -a(E- UM 2 /2 + fjN- C/D 2 / 2 ) 2 + 4T 2 (M 2 - D 2 ) , (10) 

where E = 7i is a given value of the energy, which reproduces the result achieved 
in Ref. [11]. Going to its second order version 

D= (a- 6L> 2 ) D , 

where a = AUe - AT 2 , b = 2U 2 , and e = E - UN 2 / 2 + (reduced energy), allows 
one to represent its solution in terms [12] of the elliptic cosine cn(x; k). This has the 
form D(t) = Acu[X(t - t q ); k], in which A 2 = (1 + a/X 2 )/{b\ 2 ), k = (1 + a/X 2 )/2, 
and A is a scale time that can be fixed arbitrarily [12]. 

In order to implement the requantization precedure a la Brillouin-Einstein- 
Keller |r| the complete set of the fixed point is requested. These are 

9 = 0, D = , (minimum) (11) 



6 = ±tt/2 , D = ±yjN 2 - (T/U) 2 , (maxima) (12) 
9 = ±vr/2 , D = , (saddle points if T/U < AT) (13) 
with (reduced) energy e 

e m = -TM, e M = ^(Af 2 + T 2 /U 2 ), e s = +TM , (14) 

respectively. The structure of the D — 9 phase space depends on the parameter T := 
T/MU. For r < 1/2 the domains containing the maxima (consider, for example, 
those related to 9 = tt/2 represented in Fig. la) are confined by a separatrix 



D { $ (9) = uM^l-AT 2 cos 2 {29) , (15) 

with v = +1 (u = — 1) for the upper (lower) maxima, that corresponds to the energy 
e* = Uftf 2 /2. Such disjoint domains are separated by a region for whose orbits 9 
can increase indefinitely. Another pair of separatrices confine the curves encircling 
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the minimum. These are obtained from (|8|) by setting the reduced energy e equal 
to es and read 



= ±V2N 



1 - rcos 2 (26») ± cos(2<9) WT 2 cos 2 (26l) + (1 - 2r 



1/2 



(16) 



For r = 1/2 the separatrices of the two maxima intersect in 9 = vr/2, = thus 
eliminating the possibility of ballistic evolution of the phase. At the same time (see 
Fig. lb) they merge with the minimum separatrices (|l6|). 

For r > 1/2 two-branch separatrix (pi]) forms a unique curve which is displaced 
vertically and join the upper boundary of the phase space (Dm = M) to the lower 
one (Dm = —AT)- The values in the interval 1/2 < T < 1 correspond to situations 
in which a new eight-shaped separatrix encircles the two maxima (see Fig. 2a) and 
divides the curves of motion that locally sourround each maxima from those that 
contain both. Equation ( |l6| ) still describes such a separatrix centered in the saddle 
point (Di^ = 0, 9 = vr/2), but it is convenient to introduce the variable cj) := 9 — it/2 
in view of the fact that (1 — 2r) is now negative. In the limit T — > (l/2) + the vertical 
separatrix and the eight-shaped separatrix merge, while for r — » 1 _ the maxima 
merge (see Fig. 2b) thus generating a configuration with a unique maximum which 
characterizes T > 1. 

An approximate evaluation of the energy level distribution can be easily obtained 
for curves of motion close to the extremal points just introduced. Upon expressing 
D as a function of 9 



D±(9) = vM^ 



T 2 c 2 (9) ± Tc(9)\ T 2 c 2 (9) + 1 



(17) 



UN 2 v ' W V w UN 2 
with v = ±1 and c(9) = cos(2#), the requantization procedure is enacted by setting 

I Dd9 = 2vr(n + 1/2), (18) 

where n S N, 7 is some isoenergetic closed path, and the constant fi is missing since 
D is a dimensionless action variable. 

In general, the exact calculation of the above integral is quite difficult. Therefore 
it is convenient to reduce H to a diagonal quadratic form in proximity of both 
minima and maxima. In the case of the minimum (D/J\f <C 1, 9 <C 1) the resulting 
Hamiltonians reads 

H m ~ ^- - M(T + y) + - (1 + T)D 2 + 2TM9 2 , (19) 

while the maxima entail a quadratic Hamiltonian 

Hm * UN 2 - MA+ + |[i - r~V - 2(T 2 /£y> 2 , (20) 
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24 4224 42 

e e 

Figure 1: (a) left figure illustrates the minimum separatrix (dotted lines) and the two 
separatrices (dashed lines) associated to each maximum pair for T < 1/2; (b) right figure 
shows the merging of separatrices for T = 1/2. 
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with p := [(D/M) ± (1 — T 2 ) 1 / 2 ] <C 1, <j> — 6 ± ir/2. Under such approximations 
the integral (|l^) can be easily calculated and provides the approximate form of the 
spectra 

2 



and 



E m {n) = UN 2 
E M {n) = UN 2 



UN 



r + jr(i + r) (2n + i) 



(21) 
(22) 



for the minimun and the maxima, respectively. 

The case when the Hamiltonian exhibits two maxima is quite interesting in that 
the presence of pairs of isoenergetic trajectories (spatially disjoint) suggests the 
possibility for the system to undergo a tunneling effect when considered quantum- 
mechanically. Indeed this is the case for eg < e < em, and 1/2 < T < 1 where 
isoenergetic trajectories appear around the points 



P R := (+AfVl - T 2 , ±tt/2) , P L := (-A/Vl - T 2 , ±tt/2) , 

in the phase plane (D,9). Due to the tunneling effect, close to the value T = 1/2 
the first order corrections in fi must be accounted when approximating the energy 
levels. 

A simple way to evaluate such an effect is obtained via the following standard 
scheme []l4| showing how the tunneling probability is nonvanishing when the barrier 
that separates the motions around Pr and Pl is finite. In this case a new relevant 
contribution (to the second order in the perturbative approximation) must be added 
to the energy levels (p2])- 

Let Hr and Hl be the quadratic approximation of H around Pr and Pl, re- 
spectively, and \^r) and the eigenstates of Hr and Hl corresponding to the 
same eigenvalue Em(u)\ the two eigenvalues: 



E±(n) = Em in) T # 



R 



dD 



c 



(23) 



where C is the crossing point (D,6) = (0, 7r/2), correspond then to the symmet- 
ric/antisymmetric states: |^±) = (\^r) ± \^l))/V^- At the second order we have 



$>R 



d^ R 
dD 



c 



2ttK 



cxp 



(24) 



where u = {UV^- 2 - 1)/(2T), p(D) = [2uj (hn/2 - 2ujD 2 /U)] 1/2 , and the inversion 
points D = 
in D+(0). 



±a of the motion along D have been issued from (J17|) by setting 6 = tt/2 
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4 The spin picture of the quantum dimer 



The Hamiltonian of a lattice formed by two sites (dimer) reads 
H = U(n\ + n\) + U)N - T( ai a% + a 2 af) , 



(25) 



where now N = n\ + n 2 , whereas the factor T of the hopping term in place of 
T/2 is the only difference distinguishing the open chain case from the case with 
the periodic boundary conditions a-M+j = ou with M = 2. Upon introducing the 
two-boson realization of the su(2) generators 

Jl = ^(«i«2~ + a 2af) , J 2 = ^( a i°2~ ~ a 2 af) , J 3 = i(n 2 - rii) (26) 
which satisfy the commutators 

[Jl, Jz] = iJi , [J2, J3] = iJi , [J3, Ji] = 1J2 , 

and the Casimir equation C = J\ + J| + Jf = J^J^ + 1) with J 4 = \{n 2 + ri\), 
Hamiltonian d23) can be rewritten in the nonlinear form 



H = 2[UJl - (ja + U) J± + UJl - TJ X 



(27) 



This fact has two consequences: First, the size of the Hilbert space of the system, 
encoded in the eigenvalue J of J4, is related explicitly to a macroscopic feature of the 
system since 2J4 coincides with the total particle number operator N. The latter 
naturally plays the role of constant of motion within our group-theoretic picture 
in that J4 is a group central element that is [J4, J a ] = 0, a = 1,2,3. Moreover, 
( |27| ) shows how, due the nonlinear term jf, H cannot be diagonalized via a simple 
unitary trasformation of the group SU(2). In view of the nonlinearity of ( p7[ ) the 
problem of diagonalizing H for any value of J can be faced either by expressing j| 
(and Ji) linearly in terms of the generators of a larger (dynamical) algebra, or by 
constructing approximation procedures capable of retaining the relevant features of 
the spectrum. The latter point of view is assumed in this Section. 

Consider first the low part of the spectrum of H, the features of which should 
be relevant to the zero temperature phase diagram of the BH model. The spin 
formulation of the problem allows one to implement an approximate approach based 
on the Casimir operator. In fact, the observation that the classical energy minimum 
is reached for J3 = J 2 = 0, and J\ = +VC, where C = J\ + Jf + Jf is a constant, 
leads to approximate J\ as J\ ~ VC[1 — (l/2C)(Jf + Jf)] which implies 



H 



UJl 



(p + u)Ji + ujI - tVc 



1 



1 

2C 



(Ji + Jl 



(28) 
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The variables J3, J 2 (considered "small") account for small displacements around 
the energy minimum. At the quantum level, the inequalities (J s ) <C \[C concerning 
the expectation values of J s , with C = J (J + 1) = N(N + 2)/4 >> 1, and s = 2, 3 
make natural to adopt the perturbative scheme where the operator J 2 and J3 can 
be treated as canonically conjugate variables since 

[J 2 ,J 3 ] ~ +iVC, (29) 

whereas 

[J 3) Ji] = +i^C( J 2 /VC) ~ , [J 2 , Ji] = -iVC(J 3 /VC) ~ , (30) 

are considered vanishing when compared with (29) since (J S )/VC « 1. Then, 



after labelling by the integer p the eigenstates of the harmonic oscillator term 
[1 + T/(2UVC)][Jl + W 2 J$] contained in (H), where VF 2 = [1 + 2UVC/T]' 1 , 
its eigenvalues are found to be VCW(2p + 1) so that the spectrum of H reads 



E p = ^N 2 - (/x + U)N - 2TVC + T 



(2p+l). (31) 



The approximation just implemented is of course exact for UN = 0. Hence its 
natural range of validity seems to be given by the condition T/NU > 1, which 
can be interpreted semiclassically as the requirement that the expectation value 
of the hopping term T|(Ji)| — TN/2 is greater than the (maximum value of the) 
Coulomb expectation value U(J 2 ) ~ UN 2 /4. This appears to be consistent with 
the fact that the harmonic oscillator form of the above spectrum is originated by 
the spectrum of J\ which is both linear and bounded from below. When T/NU < 1 
the spin structure of the problem fully emerges since Jf introduce the influence of 
its quadratic degenerate spectrum. 

The importance of the condition T/NU > 1 is even more evident when de- 
scribing the spectrum close to the maximum energy configuration. In this case the 
approximation J\ ~ —yC[l — (1/2C)(J| + J 2 )] makes appear in H the harmonic 
oscillator term ~[T/(2U^C) - 1][J| + 2 J|] where Q 2 = T/(T - 2U^C) which 
involves a spectrum of the form 



E q = ^N 2 - (/x + U)N + 2TVC - T 



(2q + 1) . (32) 



For T/NU < 1 the system no loger exhibits a unique maximum so that the oscillator 
approximation of the problem does not work: the squared frequency Q 2 becomes 
imaginary. The spin picture thus appears as the correct approach for understanding 
quantum mechanically the degeneracy caused by the presence of two (classical) 
symmetric maxima for T/NU < 1. 



11 



5 Two-well system dynamical algebra 



5.1 Single-boson picture 

Considering the subalgebra sujv(2), in the spin representation J = N/2 where N is 
the total number of boson, within the algebra A = su(M) allows us to recast the 
nonlinear Hamiltonian H in terms of a linear combination of generators belonging 
to A. The latter is therefore a dynamical algebra for the system. The simplest 
possible choice for A is provided by the two-boson operator realization of suq(M) 
the generators of which have the standard form 

En := b+b 3 , (i ± j) , Ha := (bfb t - b+bj)/2 (33) 

with 1 < i, j < M, and satisfy the commutation relations 

[Eij,H k {\ = - (5j k E ik -5jiEu+5iiEij-5 ik E k j) , [Eij,E kl ] = 5j k Eu-5uE k j . 

The index Q eigenvalue of the invariant operator Nb = J2i bfh ([Nb,Eij] = 0) that 
selects the specific representation of su(M) considered, denotes the total number 
of particle associated with such a bosonic realization of su(M). In fact, the dimen- 
sion of the Hilbert space basis B(M,Q) = {|ni, tim), Y^iLi n i = Q} is given by 
dimB(M, Q) = (Q + M - 1)!/[(M - 1)!Q!], where the basis states are defined as 
\ni, um) = ®iL\\ n i) an d the number operator states \rij) fulfil the equations 
bi\rii) = yjni\ni - 1) and bf\iii) = y/m + l|nj + 1). 

Setting M=N+1 (recall that N is the boson number of the two-well system) and 
Q=l one selects the fundamental realization of su(M), whose one-boson quantum 
states |0, ...,1,0...) (dim B(N + 1,1) = N+l) are in one-to-one correspondence with 
the column vector standard basis whose matrix form reads {|1) := (1,0...) T , |2) := 
(0, l...) T , ...\q) := (0, 1, 0...) T }. The generators of the resulting matrix realization 
of sui(N+l) have the form 

ll-S'jj'llgp = fiqi^jp j H-^ijIlgp = ^(^qi^ip ~ ^qj^jp) > 

where 1 < q,p < M = N + 1. The states \ J;m) constituting the standard basis of 
suat(2) can be similarly realized within the column vector basis {\q),l < q < N + 1} 
via the identification \J;m) = \q) with m = J + 1 — q. The ensuing matrix version 
of the suat(2) generators reads 

H^+llgp = \J(P~ l )( 2J + 2-p)0~ qp -l , WMlqp = (J ~P + l)Sqp 
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thus entailing that J+ and J3 can be rewritten within su(M) (recall that M = N + l) 
as 



N+l 



J + = E \M 2J + 1 

q=l 



N+l 

)E q g+l , J3 = ^2 
9=1 



J+l 



AT 



M 

Q - 2 H jq 
3=1 



(34) 



respectively. Formally, this turns out to be equivalent to representing our quantum 
description in the fully bosonic form, with 



N+l 



N+l 



J + = E v 9( 2 J + 1 - i)K h i+i > J 3 = E ( J + 1 - ?K > 

g=l q=l 

with the constraint ^2^jti bqb q = 1. Hamiltonian ( ]27| ) then becomes 

iV+l JV+l 

# = c( j) + 2c/ J2 ™ 2 (i) * q - T E J ) (^+1 + K+M > 

g=l 9 =1 



(35) 



where J) := [(J + 1/2) 2 - (771(9) - 1/2) 2 ] 1 / 2 , m{q) := J + l-q, and C( J) 
2[UJ 2 — {fx + U) J]. The related dynamics is now described by states 



N+l 

l*> = E^i°>->°> 
9=1 



(36) 



that, due to the normalization condition, turn out to define the manifold J2^=i^ l£g| 2 
in the bosonic Hilbert space. Such a manifold is an hypersphere in C N+1 and repre- 
sents the effective space in which the quantum dynamics of H can be reformulated 
in a classical form in such a way to mantain a complete equivalence with the original 
quantum problem. In fact the classical dynamics thus obtained actually concerns 
the time evolution of the components of \^f) once the initial condition |^(0)) at 
t = has been assigned. 

Based on Eq. (|36|) , the Schrodinger problem (idt — H)\^) = can be rewritten 
as a set of equations of motion 



2Um 2 ^m-T [pj (m + 1) # m+1 + Pj (m) ^ m -i] 



(37) 



where we have adopted the notation fy m := £ 9 , and p r (m) := R(q,J) with m := 
J+l — q, in order to make the present formalism closer to the spin picture introduced 
previously. Equations (37) can be derived in an independent way from the effective 
Hamiltonian 



N/2 



H[V]=C(J)+ £ [2Um 2 \^ m \ 2 -T Pj (m)(^ m -i + ^ m ^ m ) 



(38) 



m =-N/2 
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representing the energy expectation value {^>\H\^) = TC[^f], provided the Poisson 
structure {^ q ,^*} = 8 qp /ih is assumed. The eigenfrequencies associated to (38) 



and the related periodic orbits coincide with the eigenvalues and the eigenstates of 
(|37|), respectively. The latter are obtained from the solution of the secular equation 

EX m = 2Um 2 X m — T [ Pj (m + 1) X m+1 + Pj (m) X m -i] , (39) 

that is easily derived from the equations of motion ([37]). The eigenstates components 
X m in Eqs. (p9|) can be shown to be real. 



5.2 Spectrum structure 

Two symmetries characterize the structure of the eigenstates. The first is realized 
via the action of the unitary transformation U\ := exp[z7rJi] which takes J3 into 
UiJ^U^ = — J3 so that [H, Ui] = 0. The action of U\ on the standard basis is given 
by Ui\m) = exp[i7rJ]| — m). This corresponds to implementing the component 
transformation X m — > sX_ m with s = ±1 in Eqs. (|39| ) that remain unchanged 
in that Pj(m + 1) = Pj(— m) for each m. Each eigenstate has therefore a definite 
symmetry character which we make explicit by writing the eigenstates as 

\E) + = J2Sm(E)\m) 1 S m :=+S„ m , (40) 

m 

\E)_ = J2A m (E)\m), A m :=-A_ m , (41) 

m 

in the symmetric and antisymmetric case, respectively. The energy spectrum is not 
degenerate in that one has that E 7^ E' for any pair of states \E) + and \E')-. 

A further symmetry is obtained by combining the action C/3 J\ U3 = — J\ of Us := 
exp[-i7rJ3] on —TJ\ in H, with the change T — > — T which restores the initial form 
of H. The corresponding change of the eigenstate components X m — > exp[z7rm]X m 
is easily derived from the fact that the standard basis is acted on by U3 according to 
the formula U^\m) = exp[i7rm]|m). This does not make change Eqs. ( |39| ) provided 
the sign of T is reversed as well. 

An important feature of the eigenvalues structure is recognized by comparing 
the secular equations for \E) + and We first consider the case when J is a half 

of an odd integer. Then equations (|39|) can be written in the reduced form 



= (2Um 2 -E)C m -T [ Pj (m + 1) C m+1 + Pj (m) C m _ x ] , (42) 
with C = A,S, for 1/2 < m < J, while for m = 1/2 one has 

= [2U (1/2) 2 - E } C 1/2 - T \ Pj (3/2) C 3/2 + V Pj (1/2) CL 1/2 1 , (43) 
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where r] = —1 (+1) must be taken in the antisymmetric (symmetric) case. The 
remarkable consequence of the presence of r\ in (|4^) is that one may distinguish the 
eigenvalues associated with the symmetric states from those related to the antisym- 
metric ones depending on the sign of rj. How this happens is clearly shown by the 
eigenvalues equation. The latter is obtained iteratively via the recurrence equation 
for the determinant of the associated secular equation 

d[E, m] = (2Um 2 - E) d[E, m + 1] - T 2 p 2 (m) d[E, m + 2] , (44) 

which terminates with d[E, J] = 2U J 2 — E, starting from 

0=(U/2-E + V T Pj (l/2))d[E, 3/2] - T 2 p 2 (1/2) d[E, 5/2]. (45) 



Developing (|45|) by means of ( |44[) a (J + l/2)-th degree equation for E is obtained. 
For rj = no difference distinguishes the symmetric case from the antisymmetric 
case that therefore exhibit pairs of eigenstates with the same eigenvalues. As soon 
as T) is switched on each eigenvalues bifurcates relatively to the two possibilities 
r] > and 7] < 0. For small values of T the leading order in ( |45| ) is T, so that it 
reduces to 

= (U/2 -E+ V T Pj (1/2)) ® J m>1/2 (2Um 2 - E) . 

There results that only the two lowest eigenvalues are distinguished from those of the 
unperturbed case. Increasing T allows one to improve the resolution: for example, 
to the order T 2 , the four lowest eigenvalues are separated. 

The same scheme holds when J is integer (m = 0,1,2...) but the hierarchy 
described above exhibits two special cases, for m = 0, 1 : 

= -EC -Tap J (l)C 1 , (46) 

= (2U — E)C\ —T [ Pj (2) C 2 + a Pj (0) C ] , (47) 



with C = A,S, whereas Eqs. (|42| ) account for 2 < m < J. The parameter a must 
be set equal to one in the symmetric case (C = S), while in the antisymmetric case 
(C = A) the expected elimination of the component Aq is realized by setting a = 0. 
It follows that the Hilbert space dimension decreases while the secular equation will 
have a degree diminished of one. Explicitly, one has 

= —E d[E, 1] - a 2 p] (1)T 2 d[E, 2} , (48) 

which reduces to = d[E, 1] for a — > 0. In the symmetric case a ( J + l)-th degree 
equation for E issues from ( |48| ) through the same formula (|44|). 

For T sufficiently small the energy levels appear to be structured in doublets (see 
Fig. 3) that mimick the emergence of the energy degeneracy actually reached only 
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for T = 0. The eigenvalues of the doublet sector satisfy the inequality E > TN, 
consistently with the role played classically by the saddle point energy e s of formula 
(|14|). Thus the doublet number decreases when T is increased (see Figs. 4, 5). 
This effects issues from the fact that both r\ and a are multiplied by T in (^5|) and 
(45), respectively. The limit T — > thus incorporates the limits rj, a — > whereby 
the degeneracy crops up (each E is associated to a pair \E)±). The eigenstates 
generated (in the same limit) \E)± = (\m) ± | — m))/\/2 have the same self-energy 
E(0) = C(J) -\-2Um 2 of the noniteracting model (competition among the m modes 
switched off by T = 0, that is, no interwell dynamics). 

The splitting of energy levels caused by T ^ guarantees the important feature 
that eigenstates are structured so as to have components that are either symmetric 
or anti-symmetric with respect to the inversion m — > —m. A full degeneracy, in fact, 
should allow for the occurrence of states inducing a strong localization in proximity 
of those values of (J3) that classically correspond to the energy peak positions. 
In this case the eigenstates should exhibit a marked semiclassical character that 
corresponds to states where either the components with m > or those with m < 
are strongly depressed in order to permit the localization effect. 

Despite the absence of states involving a stable localization around one of the 
two energy peaks for energy values E > T/U when T/UN < 1, the realization of 
localized states can be attained by superpositions of the two states corresponding 
to the doublets. Each such a pair of states is composed by a symmetric eigenstate 
and an antisymmetric one with an energy gap AE which vanishes for T — > 0. 
The configurations where the system is localized undergo a periodic revival (with 
a life time proportional to 1/AE) which ensues from the time dependence of the 
eigenstate interference mechanism. 
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Figure 3: Energy level distribution for T — 4, U — 1, N — 41. The index m G [1,42] 
labels the energy eigenvalues. 
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Figure 4: Energy level distribution for T — 24, U — 1, N — 41. The eigenvalues are 
labelled by the index m £ [1,42]. Three doublets are still visible for m < 6 
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41: the influence of C/J| in Hamiltonian (27) is almost 
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6 Discussion 



The structure of Hamiltonian ( |38|) allows us to interpret the time evolution of a 
generic state |^} as the dynamics of an ensemble of linear oscillators with a space- 
dependent coupling (pj{m) depends on m) in the linear lattice of m modes. A 
certain amount of information can be extracted from ( |38| ) as to the form of the 
states \^>) close to the extremal configurations. 

Let us consider, first, the case T ~ 1 (we identify the classical parameter N 
with the eigenvalue N = 2J so that T = T/UN). The greatest coefficients of the 
hopping term (max[T / o J (m)] ~ UN 2 /2) are comparable with the greatest coefficient 
of the Coulomb term (max[2C/m 2 ] = UN 2 /2). On the other hand, such two classes 
of terms involve different components that are labelled by m ~ and m ~ ±iV/2, 
respectively. The minimum energy state is therefore characterized by components 
ty m = R m exp(ia m ) such that a m = a m+ i (the phase-locking involves the most 
negative the hopping term), whereas the modes \m\ <C N/2 (\m\ — N/2) must be 
maximally occupied (strongly depleted). The dynamics of wealky excited states is 
expected to exhibit a similar mode partecipation with the breaking of the (collective 
order of the) phase-locking due to small phases' oscillations. 

For opposite reasons the maximum energy state is characterized by an antifer- 
romagnetic type of phase order (a m +i = a m + n), while both the largest modes 
(|m| — N/2) and the smallest ones (|m| ^ N/2) contribute to maximize the energy. 
Such a phase order enables us to recognize the largest components of the maxi- 
mum energy state starting from Hamiltonian (|38|). Upon setting a m ~ 7r + a m _i 
to maximize the energy, and assuming that \\^ m \ — l^m-iH <C \^ m \, l^m-il which 
entails *^^ m _i + ^ m ^'^_ 1 ~ -[|^ m | 2 + |^ m „i| 2 ], Hamiltonian (||) reduces to the 
diagonal form 

N/2 

Hm~C(J)+ ]T [2Um 2 +T( Pj (m) + Pj (m + 1)) 

m=-N/2 



|\Tj I 2 



For large values of J = N/2 the label m can be regarded as a continuous parameter 
so that the highest weight coefficient [2Um 2 + T(pj(m) + pj(m + 1))] is singled 
out by differentiation. This provides m[m ± (N 2 — T 2 /U 2 ) 1 ^ 2 /2] ~ where m = 
±{N 2 -T 2 /U 2 ) 1 l 2 /2 are the values for which Y& m \ must be maximized in order to 
reach to the highest energy configuration. This result matches exactly the classical 
representation of the maxima in the phase space (see formula ([12])), except that 
quantum mechanically both the maxima are involved in the dynamics. If Y > 1, the 
solution m = must be considered which reflects in a consistent way the coalescence 
of the two (classical) maxima in a unique one. 

The regime T ^ 1 is also intersting. This case is quite sensitive to the parity 
of iV since the largest Coulomb term has gained a factor N with respect to the 
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largest T-dependent term. If N is even (i.e. m assumes integer values) the system 
ground-state must be constructed in such a way that |^o| ^> \^m\ since |^o| does 
not contribute to the Coulomb term. We thus approximate \*5?) by neglecting any 
ty m such that m ^ 0, ±1. The resulting approximation of TC[^] reads 

H[V] ~ C{J) + 4U\^i\ 2 -4T Pj (0) |*i||*o|, 

where we have assumed |^+i| = |^-i|, and the constraint 2|\& + i| 2 + |^'o| 2 = 1 holds. 
The latter suggests to set y^l^+il = sina and |^o| = cosa in the above formula, 
which in turn appears to be minimized by tga = V2T Pj (0)/U ~ TN/V2U. The 
ground-state energy thus reads 

E gs = C(J) + j[l- v/l + tg 2 a], 

that reproduces the semiclassical value E ~ C(J) — TN for T/U < 1/N. If N is 
odd (to = ±1/2, ±3/2, ...) the Coulomb term never vanishes, and the main contri- 
butions to ground-state come from *_i/2 and ^+1/2, so that the related energy is 
approximated by 

Um ^ C{J) + ^(|*_i/ 2 | 2 ± |^ + i/ 2 | 2 ) - 2T Pj (-l/2)|tt_ 1/2 ||tt +1/2 | . 

In view of the normalization l^-i^l 2 + 1^+1/2 1 = ^ reduces to 

U TN 1 

E gs ~C(J) + --—^l + 2/N, 

under the assumption |\P_i/2| = | v I / +i/ 2 |. 

The characters distinguishing the even case and the odd case suggest the inter- 
pretation of such situations in terms of insulator and supernuid states, respectively. 
In the odd case one has m = ±1/2, ±3/2, ... so that the ground-state is the super- 
position of the states \N/2, N/2 + 1) and \N/2 + 1, N/2) which involves a permanent 
one-boson tunneling between the two wells (supernuid regime). In the even case the 
ground-state essentially coincides with \ J,0) = \N/2, N/2) which entails a station- 
ary regime where the two bosonic wells are equally filled and the boson tunnelling 
is not favoured (insulator regime). It is worth noting that any contribution of the 
hopping term turns out to be definitely smaller than the Coulomb terms namely 
max[T / o J (m)] ~ U/2 < minfC/ro 2 ] ~ U when one sets F = 1/N 2 . The ensuing 
almost null importance of the hopping term is precisely the signal of the transition 
from the supernuid to the insulator regime which is known (see [|l|, [Q]), for a given 
N, to take place at T = 1/N 2 . 

A simple evaluation of the components X m of the eigenstates associated with the 
energy spectrum extremes (where relative sign of X m , X m+ \ is definite) is obtained 
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by regarding m as a continuous parameter in Equation (^). In fact, upon setting 
X m = Y m / y/(J + m)!( J — m)\ and replacing Y m ±\ with Y m ± dY m /dm in (|39|), one 
finds 

dY 

= (E - 2Um 2 + 2JT) Y m - 2Tam—^- , (49) 

dm 

where a = +1 (—1) for the minimum (maximum) energy state. The solution have 
the form Y m = Cm^exp(—Um 2 /2T) with A = E/{2To) — J and the consatant C is 
fixed by imposing the state normalization. The same treatment can be extended to 
the intermediate eigenstates but the sign oscillation of X m 's, known from numerical 
calculation of eigenstate, requires a more accurate reformulation of Eqs. ( j39|) also 
involving a second derivative term. 



7 Conclusions 

In this paper we have focused our attention on the quantum aspects of the dynamics 
of a two-site Bose-Hubbard model. The latter mimicks the interaction of two BEC's 
via tunneling effect. This has led us to study thoroughly the nontrivial structure of 
the energy spectrum the upper part of which consists of a series of doublets if T < 1, 
and has suggested the quantal interpretation of the degeneracy characterizing the 
classical configurations (a pair of symmetric orbits is associated to the same value 
of the energy) . 

The general framework in which we relate the quantum dynamics of BEC's to 
the classical formulation usually employed for macroscopic condensates has been 
described in Sec. 2, where the purely quantum model of many lattice sites (bosonic 
wells coupled within the BHM picture) is reduced to a classical lattice field theory 
by means of the TDVP procedure and the coherent state representation of the 
macroscopic (quantum) state of the lattice gas. Conversely, such a construction 
shows that a generic array of coupled BEC is naturally reconducted to the standard 



BHM and its generalization to nonperiodic lattice structures [15 when considered 
in a nonclassical regime (small or mesoscopic number of bosons per condensate). 

Sec. 3 has been devoted to illustrate the classical form of the two-site model, 
its phase space, and the orbit structure therein. After discussing the integrable 
character of the system the dynamics of which is completely accounted by canonical 
variables 6, D (see equations @), the requantization procedure a la BEK has been 
implemented. This provides integral (|l^) whose approximate solutions give the 



evaluation of energy levels plft and (p2[) in proximity of the minimum and the 



maxima, respectively. The nice feature that (22) holds for T < 1 has allowed us to 



incorporate the exchange effects via the procedure leading to formula (23). Such a 



result anticipates the basic trait of the spectrum represented by its doublet structure. 
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The reformulation of the two-site model in terms of a macroscopic spin has been 
performed in Sec. 4. Such a picture offers the possibility to approximate the energy 
spectrum in an independent way for levels close to the minimum and maximum 
values. The approximation well fits the exact energy values for T > 1 and shows 
implicitly how the exhaustive investigation of the spectrum unavoidably requires 
the spin picture. This paves the way for the construction of the dynamical algebra 
which is the base for the exact diagonalization of the model. 

The realization of s«at(2) within the larger algebra suq(N + 1) in terms of 
a single-boson representation (Q=l) has been developed in Sec. 5 and supplies 
the spin operators in the form (|34|), It provides the algebraic framework where the 



nonlinearity of J| is eliminated and the Hamiltonian becomes a linear superposition 
(see (|35l)) of sui(N + 1) generators. The procedure is easily generalized to many- 
site lattice [16|. The set-up developed leads to reformulate the two-site model as a 
(N+l)-lattice model described by the quadratic Hamiltonian (35) and a population 
constituted by a single particle. The spin nonlinear dynamics has been finally 
reconducted to a classical problem with the quadratic Hamiltonian (^) depending 
on the canonical variables (these are the components of the wavefunction solving 
the Schrodinger equation) and The proper modes of (^) coincides with the 
eigenvalues E of H that are grouped in two sectors E > TN and E < TN if 
< r < 1. The first sector contains the doublets. 

The basic result illustrated in Sec. 5 is that the levels of each doublet are associ- 
ated to a symmetric state and an antisymmetric state exhibiting almost coincident 
secular equations. Their difference, that turns out to be concentrated in one (or 
two) of the equations forming the hierarchy of equations for the eigenstate compo- 
nents (see Eqs. (|||) and jig)- (47)), has been represented via a parameter which 
joins analytically the symmetric and antisymmetric eigenvalues. This explains the 
origin of the doublet structure through a bifurcation mechanism which represents 
the main result of the present paper. In Sec. 6 the two site model ground-state has 
been interpreted for T <C 1/2 as a state with a superfluid (insulator) behavior for 
an odd (even) total particle number N thus showing how the germ of the lobe-like 
structure of the zero temperature phase diagram is already present in a two-site lat- 
tice. The method and the analysis developed in the present work can be extended 
to lattices with a number of sites greater than two. The bifurcation mechanism 
explained above is expected to provide the base for studying the chaotic character 
of the trimer case. Work on the three-site case is in progress along these lines [16|. 
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